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A GAS-KINETIC SCHEME FOR MULTIMATERIAL FLOWS 
AND ITS APPLICATION IN CHEMICAL REACTION* 

YONGSHENG LIAN* AND KUN XU* 


Abstract. This paper concerns the extension of the multicomponent gas- kinetic BGK-type scheme [26] 
to multidimensional chemical reactive flow calculations. In the kinetic model, each component satisfies its 
individual gas-kinetic BGK equation and the equilibrium states of both components are coupled in space 
and time due to the momentum and energy exchange in the course of particle collisions. At the same time, 
according to the chemical reaction rule one component can be changed into another component with the 
release of energy, where the reactant and product could have different 7. Many numerical test cases are 
included in this paper, which show the robustness and accuracy of kinetic approach in the description of 
multicomponent reactive flows. 

Key words, gas-kinetic method, multicomponent flow, detonation wave 
Subject classification. Applied Numerical Methods 

1. Introduction. The development of numerical methods for the multimaterial flows have attracted 
much attention in the past years [15, 22, 8, 9]. One of the main application for these methods is the 
chemical reactive flow calculations [2, 10, 18]. Research of reactive flows, especially the detonation waves, 
was pioneered by Zeldovich, von Neumann, and Doering, who developed a well known ZND model. The ZND 
model consists of a non- reactive shock followed by a reaction zone. Ever since the model was proposed, lots of 
theoretical and numerical work on this problem have been done. Numerical calculation of the ZND detonation 
was pioneered by Fickctt and Wood [11]. They solved the one- dimensional equations using the method of 
characteristics in conjunction with a shock fitting method. Longitudinal instability waves were accurately 
simulated. Later, Taki and Fujiwara applied van Leer’s upwind method to calculate two-dimensional traveling 
detonation waves [23, 24]. They solved the Euler equations coupled with two species equations. The chemical 
reaction was simulated by a two-step finite-rate model and the transverse instabilities around shock front 
were clearly observed. It was pointed out by Colella et al in [5] that if the numerical resolution in the 
detonative shock front is not enough, unphysical solution can be easily generated, such as the wrong shock 
speed. In order to avoid the unphysical solution, Engquist and Sjogreen [6] obtained a high order TVD/ENO 
numerical method combined with Runge-Kutta time marching scheme to solve the combustion problem and 
designed a special treatment in the shock region. Kailasanath et. al [14] extended the Flux- Corrected 
Transport (FCT) algorithm for detonations. In the early 90s, Bourlioux ct.al. combined PPM scheme with 
conservative front tracking and adaptive mesh refinement in the study of detonative waves [2, 3, 4]. They 
showed the spatial-temporal structure of unstable detonation in one and two spatial dimensions and found 
good agreement between the numerical simulation and the experimental data. Quirk [21] addressed the 
particular deficiency of the Godunov type upwind schemes in solving complex flow problems, and suggested 
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a hybrid scheme, from which he successfully simulated the galloping in one and two-dimensional detonations. 
Lindstrom [18] analyzed the poor convergence of inviscid Euler solutions in the study of detonative waves and 
suggested to solve the compressible Navier- Stokes equations directly. Most recently, Hwang et al [13] pointed 
out that not only the resolution of the reaction zone is important, but also the size of the computational 
domain is critical in the capturing of correct detonative solutions. So far, it is well recognized that a good 
scheme for the reactive flow must be able to capture correct shock speed, and resolve wave structures in 
multidimensional case, as well as present the correct period of the possible unsteady oscillation in the wave. 

Ever since the gas kinetic scheme was proposed for the compressible flow simulations [25], due to its 
robustness and accuracy it has attracted more attentions in the CFD community. In this paper, we are 
going to extend the multicomponent BGK solver [26] to high dimensions, and develop a new scheme with 
the inclusion of reactive terms. The paper is organized as follows. Section 2 introduces the governing 
equations for the chemical reactive flows in the 2D case, and describes the numerical method. Section 3 
is about the numerical experiments, which include non- reactive shock bubble interaction and ZND wave 
calculations in both ID and 2D cases. We also show a new example where the reactant and product could 
have different 7. Different from previous approach [17], the current method follows the evolution of each 
species individually, and the scheme is more robust than the previous one. 

2. Numerical Method. The focus of this section is to present a kinetic scheme to solve the following 
reacting compressible Euler equations in the 2-D case, 
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where p\ is the density of reactant, P2 is the density of product, p = p\ 4 - P2 is the total density, pE is the 
total energy which include both kinetic and thermal ones, i.e. pE — \p(U 2 + V 2 ) + Pi/fai — 1) + P2H12 — 1). 
Here U, V are the average flow velocities in the x and y directions respectively. Each component has its 
specific heat ratios 71 and 72. P = Pi 4- P 2 is the total pressure, and Qq is the amount of heat released per 
unit mass by reaction. The equation of state can be expressed as Pi — p\RT and P2 = P2RT. K(T) is the 
chemical reactive rate, which is a function of temperature. The specific form of K(T) will be given in the 
numerical example section. 

The above reactive flow equations will be solved in two steps. In the first step, the nonreactive gas 
evolution parts are solved using the multimaterial gas- kinetic method. In the second step, the source terms 
on the right hand side of Eq.(2.1) are included in the update of flow variables inside each cell. 

2.1. 2-D Multicomponent BGK Scheme. 


2.1.1. Gas-kinetic Governing Equations. The focus of this subsection is to present the multicom- 
ponent BGK scheme in two dimensions. For two dimensional problem, the governing equation for the time 
evolution of each component is the BGK model, 
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where and / ^ are particle distribution functions for component 1 and 2 gases, and g ^ and g ^ are 

the corresponding equilibrium states which and / ^ approach. The relations between the distribution 
function and the macroscopic variables are, 

J f w 4> WdS (1) + / (2 Vi 2) dH< 2 > = W 

( 2 . 3 ) = (pi,p 2 ,pU,pV,pE) T , 
where d3^ — dudvd^i, dES 2 ^ — dudvd^, 

0a * = (l,0,w,t;, i(u 2 + v 2 + ^)) T 

0 a ’ = ( 0 , 1 ,U,V, i(w 2 + l' 2 +^)) 7 ’, 

arc the moments for individual mass, total momentum, and total energy densities, £j* = Ci,i + £1,2 + * ‘ * + 
K and = £f ?1 + Cl, 2 + * • * + £2 ,k 2 * The integration elements are d£i = d£i,id£i,2-»d£i,Ki and d£ 2 = 
^£2,1 ^£2,2 • * *^£2,^2 ■ an d ^2 arc the degrees of the internal variables £1 and £2, which are related to the 

specific heat ratios 71 and 72. For the two-dimensional flow, we have 

K x -{ 5 - 371)/ (71 - 1) + 1 and K 2 = (5 - 37 2 )/(72 - 1 ) + 1 - 

Instead of individual mass, momentum and energy conservation in a single component flow, for two compo- 
nent gas mixtures the compatibility condition is 

(2.4) /(s (1) - / (1) )^ 1) <E (1 > + (g (2) - / (2) )0i 2) dH (2) - 0 , a= 1,2, 3, 4, 5. 

The equilibrium Maxwellian distributions g 11 ^ and g i2> are generally defined as 
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where Ai and A2 are function of temperature. Due to the momentum and energy exchange in particle 
collisions, in most cases, the equilibrium states and g ^ can be assumed to have the same velocity and 
temperature at any point in space and time. So, based on given initial macroscopic variables at any point 
in space and time, 

w (1) = J = {pi, P 1 U 1 , piVi, piE x ) T , 

( 2 . 5 ) W™ = J gW+WdEP' = (p 2 ,P2U 2 ,p 2 V 2 ,p 2 E 2 ) T , 

we can get the corresponding equilibrium values 
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where the common values of C7, V and A can be obtained from the conservation requirements, 
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With the definition of “averaged” value of internal degree of freedom K , 
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the values £/, V and A can be obtained from Eq.(2.7) explicitly, 
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As a result, the equilibrium states can be expressed as 
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The governing equations ( 2 . 2 ) basically correspond to viscous multimaterial governing equations, and the 
scheme presented in the next section is actually solving the Navier-Stokes flow equations, where the dissipative 
coefficients are proportional to the collision time r [25]. 

2.1.2. Multicomponent Gas-kinetic Scheme. Numerically, the Boltzmann equations ( 2 . 2 ) are solved 
using the splitting method. For example, in the x direction, we solve 


/<■’+„/<■> = £^ll, 
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and in the y direction, 
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In each fractional step, the compatibility condition (2.4) is still satisfied. 

For the BGK model, in the x direction the equivalent integral solution of / at a cell interface Xj + i /2 and 
time t is 
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for component 2, where Xi+ 1/2 is the cell interface and x f = x*_ ^ 1/2 — u(t — t r ) the particle trajectory. There 
are four unknowns in Eq.(2.15) and Eq.(2.16). Two of them are initial gas distribution functions and 
at the beginning of each time step t = 0, and the others are g ^ and g in both space and time locally 
around (x i + l/2>t - 0). 

Numerically, at the beginning of each time step t = 0, we have the macroscopic flow distributions inside 
each cell i, 

Wi = (pi , p 2 , pU, pV, pE)J . 

From the discretized initial data, we can apply the standard van Leer limiter L(5) to interpolate the conser- 
vative variables W{ and get the reconstructed initial data 
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and ( Wi(xi_i/ 2 ), ^(^+ 1 / 2 )) are the reconstructed point-wise values at the cell interfaces Xj_i/ 2 and x i+ i/ 2 . 

In order to simplify the notation, in the following x i+ i /2 = 0 is assumed. With the interpolated macro- 
scopic flow distributions W z , the initial distribution functions and / ^ in Eq.(2.15) and Eq.(2.16) arc 
constructed as 
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for component 2. The equilibrium states in Eq.(2.15) and Eq.(2.16) around (x = 0,< = 0) are assumed to be 
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and 

(2.21) 
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where H(x) is the Heaviside function, and gl'f ' 1 are the initial equilibrium states located at the cell 

interface, 
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All coefficients a\ l \a\ 2 ’ A 4 ’ ' are local constants. In order to determine all these unknowns, the BGK 
scheme is summarized as follows. 

The equilibrium Maxwellian distribution functions located on the left side of the cell interface for 

component 1 and 2 are, 
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At the location x = 0 , the relations (2.3) and (2.4) require 
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and 
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from which we have 
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Therefore, g^\ g\ 2 \ g^ and are totally determined. 

Since 77 ■ 1 ■' and <) 2 ! have the same temperature and velocity at any point in space and time, as shown 

in Eq.( 2 . 6 ), the parameters (aj 1 /' , a{^ 2 \ ajj, 2 ' , a//* ) are not totally independent. Since a\^\ oj j 2 ^, 0} ^ 

depend only on derivatives of Uo ,Vo and Ao, common velocity and temperature in space and time require 
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This is also true among the parameters . . . , a^, a^ on th c right hand side of a cell interface. So, 

inside each cell i, we have 
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The above five equations uniquely determine the five unknowns (a[*j , a{ 2 \ a* i2 , 01,3, 01,4) and the solutions is 
the following: Define 
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The solutions of Eq.(2.24) are 
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With the same method, all terms in terms can be obtained. 

By taking the limits of ( t — ► 0) in Eq.(2.15) and Eq.(2.16), applying the compatibility condition at 
(x = x i+ i/ 2 ,t = 0), and using Eq.(2.18, 2.19), we get 
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The right hand side of the above equation can be evaluated explicitly using gl^ in Eq.(2.23). Therefore, 
Pi, 0 j P 2 ,o> Ao, &o and Vo in Eq.(2.22) can be obtained from Eq.(2.25). As a result, and g ^ are totally 
determined. Then, connecting the macroscopic variables 
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at the cell interface with the cell centered values in Eq.(2.17) on both sides, we obtain the slopes for the 
macroscopic variables, 
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from which aj 1 ^ and in Eq.(2.20) and and a ^ in Eq.(2.21) can be determined using the same 
techniques for solving Eq.(2.24). At this point, there are only two unknowns A^ 1 ’ 2 ) left for the time evolution 
parts of the gas distribution functions in Eq.(2.20) and Eq.(2.21). 

Substituting Eq.(2.18), Eq.(2.19) and Eq.(2.22) into the integral solutions Eq.(2.15) and Eq.(2.16), we 
get 
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(2.27) 
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In order to evaluate the unknowns A^ 1 ' 2 ^ in the above two equations, we can use the compatibility condition 
at the cell interface Xj+j/2 on the whole CFL time step A t, 

I ( 9 (1) - f^WdE^dt + (g™ - fO^dzWdt = 0, 

from which we can get 
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73 = (1 - e~ At ' T ), 

and 

74 = — r(l - e~ At/r ) + Ate~ At/T . 

The right hand side of the Eq.(2.28) is known, therefore all parameters in /l 1 ' 1 ' 2 - 1 terms can be obtained 
explicitly. 
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Finally the time-dependent numerical fluxes for component 1 and component 2 gases across a cell in- 
terface can be obtained by taking the moments of the individual gas distribution functions and / ^ in 
Eq.(2.15) and Eq.(2.16) separately, which are 
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Integrating the above time- dependent flux functions in a whole time step Af, we can get the total mass, 
momentum and energy transports for each component, from which the flow variables in each cell can be 
updated. 


2.2. Reaction Step and Flow Update. After obtaining the flux functions across a cell interface, we 
need to solve an ODE to account for the source term, i.e. W t = S. More specifically, inside each cell we 
need to solve 


(2.29) 


' (pi)t = -K(T) Pl , 
< (P2)t = K(T) Pl , 

K (pE)t = KQoPi- 


In the current study, one step forward- Euler method is used to solve the above equations. 

In summary, the update of the flow variables inside cell (z, j) from step n to n+ 1 is through the following 
formulation, 


W n+1 = W n — 1 

^ AV 


Ay 


pAt rAt \ 

J {Fi-l/2,j — Fi+ 1/2, + Ax J 1/2 “ ^i,j+l/2)^ j 


+ A tSi y j, 


where Sij is the corresponding source terms in cell (i, j), F and G are numerical fluxes across cell interfaces 
by solving the multicomponent BGK equations, and AV is the area of the cell (z, j). 

3. Numerical Examples. In this section, we are going to test the multicomponent BGK scheme for 
both nonreactive and reactive flow calculations. For the viscous calculations, the collision time r in the BGK 
scheme presented in the last section is set to be 


T = p/P, 

where // is the dynamical viscosity coefficient and P is the total local pressure. For the Euler solutions, the 
collision time in the calculation is defined as 

r = 0.05A* + ^ ~ ^ At, 

Pi +Pr 
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where At is the CFL time step, and Pi and P r are the corresponding pressure terms in the states gi and g T 
of the initial gas distribution function /o. From the above expression, we know that in the smooth region 
there are about 20 collisions inside each time step in the current inviscid calculations, and the magnitude of 
corresponding numerical diffusion is about 1/10 of that in the Kinetic Flux Vector Splitting (KFVS) scheme 
[20, 19, 25], Also, in comparison with the previous “single component” kinetic method for the reactive flows 
[17], the current approach is more robust. The detail comparison is shown in [16]. 

3.1. Nonreactive Multimaterial Flow Calculations. In this subsection, we are going to show two 
cases about the shock-bubble interactions. The main difference between these two cases is about the initial 
density difference inside the bubble, which consequently gives different flow pattern around material interface. 


CASE(l) A M s = 1.22 shock wave in air hits a Helium cylindrical bubble 

We examine the interaction of a M s = 1.22 planar shock wave, moving in the air, with a cylindrical 
bubble of Helium. Experimental data can be found in [12] and numerical solutions using adaptive mesh 
refinement has been reported in [22]. Recently, a ghost fluid method has been applied to this case too [8]. 
A schematic description of computational set-up is shown in Fig. (4.1), where reflection boundary conditions 
are used on the upper and lower boundaries. The initial flow distribution is determined from the standard 
shock relation with the given strength of the incident shock wave. The bubble is assumed to be in both 
thermal and mechanical equilibrium with the surrounding air. The non-dimensionalized initial conditions 
are, 

W = (p — 1 , U = 0, V = 0, P = 1,7 = 1.4), pre-shock air 
W = (p= 1.3764, U = -0.394, V = 0, P = 1.5698 , 7 = 1.4), post-shock air 
W = (p = 0.1358, U = Q,V = 0,P = 1 , 7 = 1.67), Helium. 


In the computation, the nondimensional cell size used is Ax = Ay = 0.25. 

In order to identify weak flow features which are often lost within contour plots, we present a number 
of Schlieren images. These pictures depict the magnitude of the gradient of the density field, 


(3.1) 


1^1 = V ( i )2+( ! )2 ’ 


and hence they may be viewed as idealized images; the darker the image the larger the gradient. The 
density derivatives are computed using straightforward central-differencing. The following nonlinear shading 
function, <f> is used to accentuate weak flow features [22], 


(3.2) 


(p = exp (—k 


|Ap| 

|Ap| mox ; ’ 


where A; is a constant which takes the value 10 for Helium and 60 for air. For R22 simulation in the next 
test case, we use 1 for heavy fluid and 80 for air. 

Fig. (4.2) shows snapshots of Schlieren-type images at nondimensional time t=0.0 and t=125.0. Before 
the shock hits the bubble, wiggles usually appear around the bubble because the numerical scheme cannot 
precisely keep the sharp material interface. The wiggles spread in all directions. When they reach the 
solid wall, they bounce back. But all these noise have a very small magnitude. After the shock hits 
the bubble, the original shock wave separates into a reflected and a transmitted shock waves. A complex 
pattern of discontinuities has formed around the top and bottom of the bubble. Since Helium has a lower 
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density in comparison with air, any small perturbation at the material interface can easily be amplified to 
form the instability. This instability at the material interface is closely related to the Richtmyer- Meshkov 
instability. In comparison with the result in [8], the current scheme could capture the unstable interface 
structure automatically. The result here is basically consistent with both the experiment and that from the 
mesh-refinement study [22]. It is an interesting problem to further study shock-bubble interaction case, and 
understand the dynamics of any special numerical treatment on the interface stability. In our calculations, 
we do not specifically pick up the location of interfaces. 

CASE(2) A M s — 1.22 shock wave hits a R22 cylindrical bubble 

With the same scheme, we investigate the interaction of a M s — 1.22 planar shock wave, moving in the 
air, with a cylindrical bubble of R22. The main difference between this case and the previous one is that the 
density of the bubble here is much larger than the density of air. The initial data is as follow 

W = (p = 1, U = 0, V = 0, P — 1, 7 = 1.4), pre-shock air 
W = {p= 1.3764, U = -0.394, V = 0, P = 1.5698, 7 = 1.4), post-shock air 
W = (p = 3.1538, U = 0, V = 0, P = 1, 7 = 1.249), R22. 

In the numerical experiment we use Ax = Ay = 0.25. Fig. (4.3) shows two snapshots of Schlieren-type images 
at nondimensional time t=0.0 and t= 150.0. Due to the higher density in the bubble region, different from 
Case (1) the material interface in this case is basically stable. This observation is also consistent with the 
theoretical understanding and physical experiment. 

3.2. Reactive Flow Calculations. The study of detonation wave has been undertaken theoretically 
and computationally for over a century. The successful theory of ZePdovich, von Neumann, and Doering has 
come to be a standard model. The ZND solution for the reacting compressing Euler equations is described in 
[10], which consists of a non- reactive shook followed by a reaction zone; both the shock and the reaction zone 
travel at a constant speed D. Given 7 and heat release Qo> there is a minimum shock speed, the so-called 
Chapman- Jougnet value Dqj , above which the ZND solution can be constructed. 

The parameter which relates to the shock speed D of a given detonation wave to the C J velocity Dc j 
is the overdrive factor f, which is defined as 



The value of f determines the stability of the detonative front. 

In the following test cases, we only consider the reactive flow with two species, i.e. the reactant and 
the product. The reactant is converted to the product by a one-step irreversible reactive rule governed by 
Arrhenius kinetics. The factor K(T), which depends on the temperature, is given by 

K(T) = K 0 T a e- E+/T , 

where Kq is a positive constant. In the current paper, we assume that a = 0 and the gas constant R is 
normalized to unity. Therefore, the above temperature T is determined by T = P/p. 

One important parameter in the numerical calculation of ZND solution is the half- reaction length L\j 2 , 
which is defined as the distance for half-completion of the reactant starting from the shock front. Usually 
the reaction prefactor K 0 is selected such that the half-reaction length is unity. From the Arrhenius formula, 
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the half reaction length is defined as 


(3.4) 



D — U 

K 0 Zgx p(^-) 


dZ, 


where D is the speed of the shock, U is the post-shock flow speed. 

In the output of numerical results, the mass fraction Z is defined as 


Pi + p2 * 

Case(l): 1 -D stable ZND detonation: 7 = 1 . 2 , Qo ~ 50, E+ = 50.0, f = 1.8 

This test case is from [ 2 ]. The pre-shock state is normalized to Po = Po = 1 and velocity Uo = Vo = 0, 
the post-shock can be obtained using Chapman-Joguet condition. The prefactor K 0 is chosen to be K 0 = 
145.68913 so that the length of the half-reaction zone L \/ 2 is unit. This case corresponds to the stable ZND 
profile. The results with 10, 20 and 40 points/L^ are shown in Fig. (4.4) and (4.5). 


Case( 2 ): 1-D unstable detonation: 7 = 1.2, Qo = 50, E + = 50, f = 1.6 

In order to get a high quality simulation result for the unstable overdriven detonation, a high resolution 
solution is usually required to resolve the instability. At the same time, the correct capturing of oscillatory 
period requires a large computational domain. As pointed out in [13], for a particular computation, one can 
be tempted to keep only a few points behind of the shock, with the reasoning that the information behind 
the shock either never catches up with, or does not affect the shock during the computation. However, if 
too small a computational domain behind the shock is specified, the points at the edge of and outside of 
the computational domain cease to be updated after some time, leading to a corruption of the data in that 
region. The U + c waves emanating from inappropriate boundary condition eventually catch up with the 
shock itself, thus erroneously alternate the shock properties. The analysis in [13] shows that if one expects 
the numerical results at time t to be correct, the computational domain L and t must satisfy the following 
inequality 


(3.5) 


t < 


L L 

U + c- D + D' 


where U is the speed of the post-shock flow, and c is the sound speed. For the current test, L should satisfy 


L > 1.88*. 


This classical unstable detonation wave was first used by Fickett and Wood [ 11 ]. An important physical 
quality for unstable detonation is the pressure history at the precursor shock in the oscillatory ZND wave 
as a function of time. For a stable ZND wave, this shock pressure history should exhibit small fluctuations 
about the known precursor shock value and decay as time evolves. In the case of unstable detonations, the 
shock front pressure history makes larger excursions from the ZND value. For the case with 7 = 1.2, qo = 50, 
E + = 50, and overdrive f = 1 . 6 , according to Erpenbeck [7] this ZND profile is a regular periodic pulsating 
detonation with maximum shock pressure per period given by 101.1 ± 0.2 while the unperturbed ZND shock 
pressure is 67.3. 

In the current study, the density and pressure are normalized to unit after shock. According to Qo = 
50 , 7 = 1.2, the CJ speed becomes Dqj — 6.80947, and the prefactor is chosen to be Ko = 230.75 so as to get 
unit half-reaction length. The post-shock state can be determined by Chapman- Jouguet condition with the 
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Table 3.1 


Maxima and minima pressure vs. time for f= 1.6 case and 80 /Lij 2 - 


Time 

Maxima 

Time 

Minimum 

7.3513 

114.1553 

11.8038 

60.1576 

15.9353 

85.0627 

18.9221 

56.7383 

23.3201 

98.1318 

26.3057 

56.7478 

30.7833 

98.3344 

33.6993 

56.8976 

38.1373 

97.8645 

41.1103 

56.7854 

45.6102 

98.0387 

48.6158 

56.5972 

53.1075 

98.8378 

56.0587 

56.8738 

60.5059 

98.1242 

63.4607 

56.9737 

67.9318 

97.3600 

70.8918 

56.6064 

75.4233 

98.6184 

78.3885 

56.6841 

82.8773 

98.7023 

85.8014 

57.0227 

90.2201 

97.3901 

93.2212 

56.7298 

97.6928 

98.2211 




given shock speed. Due to the “start-up” numerical incompatibility, there is a large initial shock pressure 
up to 114 at time t equal to 8, see Fig(4.6). After t > 15, the motion of the shock front becomes periodic. 

In this test, we observe that at least 20 points/Z^/2 is needed for getting a correct unstable ZND solution. 
In Fig. (4.6) and (4.7) we show the numerical results with 20 points/Li/2 and 40 points/Li/2 respectively. 
At the same time, the result with 80 points/Li/ 2 is given as a reference. In Table(3.1), the data of local 
maximum and minimum pressure as a function of time are listed. 

Case (3) Weak shock wave hitting the reactant 

In order to validate the multicomponent BGK scheme, we design the following ID case to simulate the 
chemical reaction in which the reactant and product have different 7. The initial condition is given below, 

Wl = ( pi , Ul, Pl, 7 l) = (2.667, 1.479, 4.500, 1.4) ' post-shock air 
W M = {Pm Mm, Pm, 7m) = (1.0, 0.0, 1.0, 1.4) pre-shock air 
W R = (pR, U R , P r , 7^) = (0.287,0.0,1.0,1.2) (Reactant). 

This case is about a weak shock wave with M = 2.0 hitting the reactant. We use the Arrhenius form for the 
reaction rate with E+ = Q 0 = 50, and Kq = 600.0. The numerical cell size is Ax = 1/2000. In Fig. (4.8) we 
show the numerical results at time t = 0.20. Since the shock is too weak to construct a ZND solution, the flow 
motion looks only like a two-component nonreactive gases. From Fig. (4.8), we can see the ordinary incident 
shock moves faster than that of the transmitted shock, and the weak reflection wave moves backward. 

Case (4) Strong shock wave hitting the reactant 

We increase the strength of the shock in Case(3) up to M — 8.0. The initial condition is given below, 
Wl = (pL, Ul, Pl^l) = (5.565,7.765,74.50,1.4) post-shock air 
W M = (pmi^a/i Pm ,7m) = (1.0, 0.0, 1.0, 1.4) pre-shock air 

W R = (p R , Ur , P R , 7*) = (0.287, 0.0, 1.0, 1.2) (Reactant). 
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In Fig. (4.9) we show the numerical results at time t ~ 0.05. From the figure, we observe that after the shock 
hits the reactant, a ZND solution is obtained. 

Case(5) Viscous Reactive Flow 

This case is from [18]. The initial data is a one-dimensional ZND profile in the x-direction. The ZND 
wave connects the left state pi = 1.731379, U\ — 3.015113 Vj = 0, piEi = 130.4736 by a Chapman- Jouget 
detonation with the right state p r = 1, U r = 0, Vr = 0 ,p r E r = 15. If no transverse gradient is present 
in the initial data, the numerical scheme will preserve the one- dimensional ZND profile. So, a periodic 
perturbation is imposed in the y-direction of the initial ZND profile, where the initial data W(x,y, 0) is set 
to Wznd{x + Aa:NINT(^^cos(47ry))), where NINT(z) is the nearest integer to z. 

The current test has Qq = E + = 50, 7 = 1.2. The reaction rate Kq is set to be 10 4 . The coefficient of 
dynamical viscosity p is set to 1.0e-4. With the above choice of parameters, the half reaction length L\/ 2 
of the inviscid one-dimensional Chapman-Jouget detonation wave is equal to 0.0285. In our computation, 
Ax = Ay = is used. Therefore, there are about 23 points/Li/2- 

Based on the analysis in [13], in order to get an accurate solution it is sufficient to use a computational 
domain x e [0,1.2]. At the left and the right boundary, we impose the left and right states of the initial 
traveling wave solution. At the lower and upper boundaries, periodic boundary conditions are used. 

Fig. (4. 10) shows a sequence of snapshots of the density distributions starting from the time t = 0.0. 
Fig. (4. 11) is the snapshot of pressure at later times when the ZND front has a regular periodic oscillating 
profile. The first picture is taken at t = which is just after the collision of two triple points. This figure 
clearly shows the formation of a Mach stem. In the next few snapshots, the movement of triple points along 
the transverse shock front are clearly captured. A high pressure spot develops at the location of triple-point 
intersection. Fig. (4. 12) shows the snapshots of the temperature variations. More figures, such as the mass 
fraction and vorticity, are presented in [16]. 

This test case corresponds to the cellular regime [2]. The hot spots in the shock front should display a 
regular diamond propagating pattern, such as observed in physical experiments. In Fig. (4. 13), we plot the 
numerical soot track of the location of shock front, which is the successive geometric representation of the 
ZND front profile. Since only the position of the shock front is recorded, one dimensional data is required 
at each output time. From the numerical soot track display, we can easily observe the formation of cellular 
pattern. 

4. Conclusion. In this paper, we have successfully extended the BGK-type gas- kinetic scheme to 
multidimensional reactive flows. Since each component of the flow is captured individually, mass conservation 
is precisely preserved for each component in nonreactive multimaterial flow calculations. For the reactive 
flow calculations, the mass exchange between different components has been implemented in the current 
kinetic method, as well as the energy release in the reaction process. Many numerical test cases validate 
the current approach and show the advantages of the kinetic scheme in the description of multicomponent 
flow calculations. For example, the unstable and stable material interfaces are captured automatically in the 
shock-bubble interaction cases. 
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Fig. 4.2. Numerical Schlieren images of the interaction of an M s =1.22 shock wave in the air moving from right to left 
over an Helium cylindrical bubble. The second image is the density distribution at time t = 125.0. The third one is the density 
profile along the centerline of the second figure. 
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Fig. 4.3. Numerical Schheren images of the interaction of an M s = 1.22 shock wave moving from right to left over an R22 
cylindrical bubble. The second image is the density distribution at time t = 150.0. The thrid one shows the density profile 
along the centerline of the second figure. 
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Fig. 4.5. Numerical solutions (solid line) of density p, velocity U, pressure P and mass fraction Z, where f=1.8, 7 = 1.2, 
Qo = E + = 50, L 1/2 = 10, and 10 pomts/Li / 2 (CFL—0.5). The dash line is the exact solution. 
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Fig. 4.7. Local maximum pressure variation as a function of time for the overdriven detonation, where f=1.6, 7 = 1.2, 
Q 0 = E + = 50, and L^j 2 = 1.0. Solid line: 80 points dash-dot line: 4 0 points/L\f 2 (CFL=0.5). 
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Fig. 4.8. Weak shock wave (M = 2.0) in the air ('y = 
has E+ = Q 0 = 50, and K 0 = 600.0 (CFL=0.5). 
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Fig. 4.9. Strong shock wave (M = 8.00,) in the air (7 = 1.- 
has E+ = Qq = 50.0., and K 0 = 600.0 (CFL=0.5). 
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Fig. 4.12. Sequence of ten snapshots of temperature start from time t = || with time increment of where Qq = E + = 


50,7 = 1*2, Ax = Ay = gfer, 23 points/L 
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Shock moves from left to right. 
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Fig. 4.13. Numerical soot track: successive geometric representation of the ZND shock front. Each line represents the 
location of shock front. The vertical axis corresponding to the time. 
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